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ABSTRACT 

Direct multi-dimensional numerical simulation is the most reliable approach for calculating the fluid dynam- 
ics and observational signatures of relativistic jets in gamma-ray bursts (GRBs). We present a two-dimensional 
relativistic hydrodynamic simulation of a GRB outflow during the afterglow phase, which uses the fifth-order 
weighted essentially non-oscillatory scheme and adaptive mesh refinement. Initially, the jet has a Lorentz factor 
of 20. We have followed its evolution up to 150 years. Using the hydrodynamic data, we calculate synchrotron 
radiation based upon standard afterglow models and compare our results with previous analytic work. We find 
that the sideways expansion of a relativistic GRB jet is a very slow process and previous analytic works have 
overestimated its rate. In our computed lightcurves, a very sharp jet break is seen and the post-break lightcurves 
are steeper than analytic predictions. We find that the jet break in GRB afterglow lightcurves is mainly caused 
by the missing flux when the edge of the jet is observed. The outflow becomes nonrelativistic at the end of the 
Blandford-McKee phase. But it is still highly nonspherical, and it takes a rather long time for it to become a 
spherical Sedov-von Neumann-Taylor blast wave. We find that the late-time afterglows become increasingly 
flatter over time. But we disagree with the common notion that there is a sudden flattening in lightcurves due 
to the transition into the Sedov-von Neumann-Taylor solution. We have also found that there is a bump in 
lightcurves at very late times (~ 1000 days) due to radiation from the counter jet. We speculate that such a 
counter jet bump might have already been observed in GRB 980703. 
Subject headings: gamma-rays: bursts - hydrodynamics - methods: numerical - relativity 



1. INTRODUCTION 

In the standard fireball shock model for gamma-ray bursts 
(GRBs), afterglows are due to synchrotron emission pro- 
duced during the slowdown o f GRB outflows by surround- 
ing media (see e.g., Zhang & MeszarosI 120041 : iPiranI 120051 : 
lMeszarosll2006t lGranoti2007i for recent reviews). The out- 
flows of GRBs are believed to be ultrarelativistic jets. The 
deceleration of the ultrarelativistic jet-like outflow in the 
early afterglow stage can be well described by the spherical 
Bland ford-McKee self-similar solution (Bla ndford & McKea 
11976*). This is because the relativistic beaming effect makes 
the material in the jet behave like an angular patch of a spher- 
ical blast wave. As it sweeps up the surrounding medium, the 
jet decelerates and becomes less and less relativistic. Hence 
the jet-like outflow will eventually expand sideways and be- 
come increasingly spherical. Thus the GRB outflow in the 
very late stages can be described by the Sedov-von Neumann- 
Taylor non-relativistic self-similar solution. Unfortunately, 
there is no good analytic solution that can describe the dy- 
namics of the sideways expansion and the transition from the 
ultrarelativistic phase to nonrelativistic phase. 

The evolution of GRB outflows is an extremely important 
problem. The modeling of observable afterglow emission 
depends upon the dynamics of the outflow. An achromatic 
break in afterglow lightcurves is observed in some GRB s 
(e.g., GRB 990510. lHarrison et alJ[T999l: IStanek et al][T999h . 
The iet break is an indication that the GRB o utflow is jet-Uke 
(lRlioadslll999l: ISari. Piran. & Halpernl ll999^. To understand 
the jet break, we must understand the multi-dimensional dy- 
namics of a relativistic jet. The late-time radio afterglow can 
be useful for an accurate estimate of the total energy of the 
GRB outflow dFrail et alJ bOOO: Ber ger et allllOOllFrail et all 
l2005h . Thus, the understanding of the transition from an ultra- 
relativistic jet into a Newtonian spherical blast wave is crucial 
for correctly interpreting the late-time radio observations. 



The dynamics of GRB outflows is a multi-dimensional 
problem except during the early Blandford-McKee and late 
Sedov-von Neumann-Taylor phases. Analytic approaches 
have been attempted to model the sideways expansion based 
upon the simplified model of the jet as having a top-hat struc- 
ture as a function of angle during the expansion ( Rhoad^ 
ll999HSari et alJl999l:lPanaitescu & MeszarosI 19991) . In fliese 
models, the jet experiences significant sideways expansion. 
To include multidimens ional effects in an inexpensive way, 
iKumar & GranotI (l2003h reduced the multi-dimensional hy- 
drodynamic equations to one dimension by integrating over 
the radial profile of the jet. They found little sideways 
expansion of the jet. However, the validity of these ap- 
proaches should be tested against full multi-dimensional hy- 
drodynamic simulations, which is certainly the most reli- 
able approach. Unfortunately high-resolution multidimen- 
sional relativistic hydrodynamic simulations are very expen- 
sive. Thus far, very few multidimensional hydrodynamic sim- 
ulations have be en performed to inv estigate the dynamics of 
GRB outflows. iGranot et alJ (1200 ll) in their pioneering nu- 
merical work found that GRB jets experience very little side- 
ways expansion. But their simulation was not long enough to 
cover the transition into the Sedov-von N eumann-Taylor solu- 
tion. A recent work by lCannizzo et al.l (|2004) has also found 
very little expansion in a different initial setup, but their sim- 
ulation suffered from severely low resolution. 

In this paper, we present a high-resolution relativistic hy- 
drodynamic simulation of the evolution of a GRB outflow 
during the afterglow phase. The evolution spans from when 
the Lorentz factor of the jet is 20 until 150 years after the 
burst. This high -resolution simulation is possible because 
our code, RAM (IZhang & MacFadyenI 12006 ). is massively 
parallel, utilizes adaptive mesh refinement (AMR), and uses 
a fifth-order method. The initial setup and results of the 
hydrodynamic simulation will be presented in Section |2l 
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Based u pon the hydrodynamic data and standard afterglow 
mode ls (ISari. Piran. & Narayaiil 119981: iGranot. Piran. & Saril 
Il999l) . we have calculated synchrotron emission from the sim- 
ulated outflow (§[3]l. We conclude with further discussions of 
our results (§|4|l. 

2. DYNAMICS OF GRB OUTFLOWS 

2. 1 . Numerical Method and Initial Setup for Hydrodynamic 
Simulation 

The special relativistic hydrodynami c simulation in this pa- 
per w as performed with the RAM code (Zhws. & MacFadyenI 
120061) ■ RA M utihzes the A MR tools in the FLASH code 
version 2.3 ("Fryxell et al. 2000), which in tu rn is a modified 
versio n of the PARAMESH AMR package (iMacNeice et all 
l2000l) . RAM includes several modules to solve the spe- 
cial relativistic hydrodynamics equations. In the simu- 
lations in this paper, the fifth-order weighted essentially 
non-oscillatory (WENO) scheme was used. This scheme 
has be en shown to achieve fifth- order accuracy for smooth 
flows dZhang & MacFadyenll2006l) and excellent treatment of 
shocks and contact discontinuities with no tunable numerical 
parameters. 

It is common in hydrodynamic simulations to use a constant 
gamma-law equation of state (EOS) given by 

P = (l-l)pt, (1) 

where P is pressure, p is mass density, e is specific internal en- 
ergy density and 7 is the adiabatic index all measured in the 
fluid rest frame. The adiabatic index can be assumed to be 4/3 
and 5/3 for relativistically hot and cold gas, respectively, but 
it is usually set as a constant for entire simulations. The prob- 
lem we studied in this paper, however, involves both hot and 
cold gas. GRB outflows are relativistically hot at the begin- 
ning of the afterglow stage and become Newtonian in terms of 
both fluid speed and sound speed. The effect of the adiabatic 
index on GRB afterglows can be quite large. For example, the 
density jump due to Newtonian strong shock compression is a 
factor of 4 and 7 for 7 = 5/3 and 4/3, respectively. Hence the 
constant gamma-law equation of state should be avoided for 
accurate treatment of afterglow dynamics. The exact equation 
of state for ideaXj;as, which works for arbitrary temperature, 
was given by fSyngd (Il97l1) . However, it is very expensive 
to use in numerical simulations because it involves modified 
Bessel functions of the second and third kinds. To accurately 
follow the evolution of GRB outflows with a reasonable cost, 
we used the TM EOS proposed by Mignone, Plewa, & Bodo| 
(^05). This equation of state has the correct asymptotic be- 
havior to 7 = 5/3 and 4/3 in the nonrelativistic and relativistic 
temperature limits, differs with the exact equation of state by 
less than 4% and can be solved at negligible cost. 

Our initial model is derived from the Blandford-McKee 
solution ( Blandford & McKed Il976h . Given the energy of 
the blast wave, Ei^o and the density of the cold surrounding 
medium, po, the blast wave can be fully described by a set 
of self-similar relations that give the Lorentz factor, pressure 
and density in the shocked fluid. In our simulation, the en- 
ergy in the Blandford-McKee solution is set to E^^o = 10''^ erg. 
Note that this energy is not the total energy in the outflow. 
Instead, it is the compensated energy the outflow would have 
if it were a spherical blast wave. In the simulation in this 
paper, we adopted a number density of Icm"^ for the sur- 
rounding medi um, which is consistent wit h GRB afterglow 
modeling (e.g., iPanaitescu & Kuinaill200lL 12002) . The sur- 
rounding medium is assumed to mainly consist of protons. 



and it is assumed to be cold, Pq <^ pQ. Numerically, we set 
the pressure to Pq = lO'^'po, where units in which the speed 
of light c = 1 are used in the simulation and henceforth in this 
paper The hydrodynamic simulation is started at the moment 
when the Lorentz factor of the fluid just behind the shock is 
70 = 20. We have run a two-dimensional axisymmetric simu- 
lation, in which the half opening angle of the GRB jet is set 
to be 9q = 0.2 radians. Thus the total energy in the twin jet is 
Ej ~ 2.0 X 10^' erg. According to observations, a total energy 
of a few times 10^' erg is typical (e.g.. iFrail et al.ll200li) . A 
half opening angle of 0.2 radians is also reasonable because 
it is within the r ange of half opening angle inferred from ob- 
servations (e.g.. Frail et aLlUoOU IPanaitescu & Kumaij|2001l 
|2002; Zeh et al. 2006). We chose a slightly large opening an- 
gle because of numerical reasons. It is extremely expensive to 
start our simulations from very early times when the Lorentz 
factor of the jet is larger, because most of the energy in the 
Blandford-McKee solution is concentrated in an extremely 
thin layer behind the shock with a width of ^ R/j^ in the lab 
frame, where R is the radius of the shock. The "lab" frame is 
the frame in which the central engine is at rest and is equiv- 
alent to the observer frame ignoring cosmological factors. In 
our simulation, we need to resolve the thin structures behind 
the shock in order to accurately capture the true dynamics. 
Furthermore, we want the opening angle and Lorentz factor of 
the jet to satisfy the relation, 70 ^ 1 /9q, so that the Blandford- 
McKee solution is still valid for the jet. Our choices of 70 = 20 
and 6*0 = 0.2 are thus reasonable. 

The initial radius of the blast wave at the beginning of the 
simulation is Rq ~ 3.8 x lO'^cm. The Sedov length for an 
explosion with an energy of Ej ~ 2.0 x 10^' erg into a medium 
with a density of po = 1-67 x 10~^'*gcm"-' is 

/ E V'^ 
km = ,, ' 2 ^ 6.8 X lO'^cm. (2) 
V(47r/3)p()cV 

To simulate the evolution of the blast wave up to the New- 
tonian regime, a large numerical box is necessary. The size 
of the numerical box for our two-dimensional simulation is 
^MAX = 1.1 X 10'^ cm, which is about 16 Sedov lengths. Most 
of the energy in the blast wave is initially concentrated in 
a thin layer behind the shock with a width of A ~ 10'^ cm 
in the lab frame. It is thus very challenging to simulate the 
blast wave over such a large dynamic range in lengthscale 
even with AMR. At the beginning of the simulation, the thin 
shell behind the shock occupies a relatively small volume. 
Thus we can afford to use very high resolution initially. As 
the blast wave expands, the volume of the thin shell, which 
contains most of the energy, becomes larger. Thus AMR 
needs to create increasingly more numerical cells in order 
to maintain the same high resolution of the structure, mak- 
ing the simulation prohibitively expensive. In our simula- 
tion, an alg orithm for de refinement of the adaptive mesh sug- 
gested by GranotI (l2007h was used to save computing time. In 
the FLASH/PARAMESH AMR package, there is a parameter 
which controls the maximal level of refinement and therefore 
decides the size of the finest grid, which is usually fixed. We 
change this parameter over time according to an algorithm 
which utilizes a nice property of the blast wave that the thin 
shell which needs to be resolved during the relativistic regime 
behaves as A ^ 7?/7^ ^ f"*, where t is the time in the lab 
frame. Our algorithm reduces the maximal refinement over 
time so that roughly the same number of cells are used in the 
r-direction to resolve the radially widening shell at different 
times without suffering from low resolution. 
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A spherical grid (r,9) with < r < 1.1 x lO'^cm and < 
< 7r/2 is employed in our two-dimensional simulation. Ini- 
tially 16 levels of refinement is used, and the finest cell has 
a size of Ar ~ 5.6 x 10'^ cm and A0 ~ 9.6 x 10"\ At this 
resolution, there are 2086 cells at (^-direction inside the jet, 
and 17 cells inside a shell with a width of Aq = Rq/jq be- 
hind the shock front. If a uniform grid were used, the total 
number of cells would have been more than 3 billion in or- 
der to achieve the highest resolution provided by AMR. With 
AMR, less than 7 million cells are needed. In our derefine- 
ment algorithm, the maximal level of refinement decreases in 
response to the change of the shell width until a specified min- 
imal numerical resolution is reached. The minimal resolution 
is chosen such that AMR uses at least 1 1 levels. 

2.2. Results of Two-Dimensional Hydrodynamic Simulation 

There are several time scales that are commonly used in the 
literature. These time scales are measured in the lab frame. 
Note that they are different from observer times. 

• tg: This is the time at which the Lorentz factor behind 
the shoc k is equal to 7 = 1 /9q for a Blandford-McKee 
solution (iBlandford & McKeelll976h . It is given by 

to « 373Eli',,n-'^' (^^) ' days, (3) 

here £150,53 is the isotropic equivalent energy Ziiso in 
units of 10^^ erg and no is density of the medium in units 
of cm~^ . 

• ts'. iLivio & WaxmanI (120001) argued that the transition 
from a jet to a spherical self-similar solution takes place 
over a time Af, « Rg /c where Rg is the jet radius at 
time tg. Thus the outflow would become spherical at fj, 
which is given by 

_ / /, \ 2/3 
/, « tg +Rg/c « 745£i'/„'53«-'/' f ^ j days. (4) 

• tsm- iLivio & WaxmanI (|2000|) defined the time fsNT 
as the time at which the shock front moves at the 
speed of light assuming the blast wave is the Sedov- 
von Neumann-Taylor solution with an isotropic explo- 
sion energy of Ej. They argued that the GRB out- 
flow becomes subrelativistic and can be described by 
the Sedov-von Neumann-Taylor solution after the time 
fsNT, which is given by 

Note that the time fsNT is earlier that the initial time 
of our simulation to, when the Lorentz factor is 70 = 20 
for a relativistic jet obeying the Blandford-McKee solu- 
tion. The jet has a Lorentz factor of 29 at Jsnt- Thus, the 
GRB outflow at f = fsNi is far from being described by 
the Sedov-von Neumann-Taylor solution. The discrep- 
ancy is caused by assuming that the outflow is isotropic. 

• fNR: iPiranI (l2005 l) argued that the transition from rela- 
tivistic to Newtonian should take place at Jnr = ^nr/c, 
where /nr = (3£'iso/47rp()C^)'''^ is the Sedov length as- 
suming that the GRB jet does not expand sideways. Af- 
ter the transition the GRB outflow can be described by 



the Newtonian Sedov-von Neumann-Taylor solution. 
The time is given by 

fNR «970£//„'g3«o'/' days. (6) 

At the time Jnr, the Lorentz factor behind the shock is 
1 .2, if one assumes that the GRB does not expand side- 
ways and the Blandford-McKee solution still applies. 

It should be emphasized again that the above time scales are 
measured in the lab frame. They are different from observer 
times unless the GRB outflow is already nonrelativistic. We 
will further discuss these time scales in Section ITTI 

Figures [T] and |2] show a series of snapshots of the two- 
dimensional hydrodynamic simulation. The simulation starts 
at to ~ 147 days (measured in the lab frame of the burster) 
after the initial explosion, and it was run up to 150 years 
after the explosion. The surrounding medium has a den- 
sity of 1.67 X 10~^^gcm"^ and a specific internal energy of 
2.25 X 10~'"'ergcm"^. Initially the relativistic outflow prop- 
agates mainly along the radial direction. Later the jet will 
inevitably undergo sideways expansion. The structure of the 
flow is very complicated. In Panels (b) and (c) of Figures [T] 
and |2] we can identify a shock moving sideways, a rarefac- 
tion wave propagating towards the jet axis, and a contact dis- 
continuity in between. The contact discontinuity is Kelvin- 
Helmholtz unstable. There is also a reverse shock inside the 
rarefaction wave propagating towards the axis. This is caused 
by the deceleration of the material which moves sideways and 
sweeps up more and more surrounding medium. It is shown 
in Panels (b) and (c) of Figures [T] and |2] that the morphology 
of the GRB outflow at early times consists of a jet cone which 
moves primarily in the radial direction and a surrounding lobe 
which moves both radially and sideways. At late times, the 
GRB outflow becomes egg-like and then increasingly spher- 
ical (Panels (d), (e) and (f) of Figures [T] and |2]l. The nearly 
vertical Mach stem at the equator is due to the collision of 
shocks from opposite hemispheres. This feature is unlikely to 
have direct observational consequences, but it is an indication 
of the accuracy of the simulation. 

What is important for observations is the distribution of 
the bulk of the energy. Since the jet material near the for- 
ward shock front initially moves at nearly the speed of light, 
very little sideways expansion takes place for the ultrarela- 
tivistic material due to relativistic kinematics. Inside the jet, 
the Lorentz factor decreases radially inwards. The mildly rel- 
ativistic and Newtonian jet material behind the forward shock 
front undergoes more sideways expansion as shown in Fig- 
ures [T] and |2] Also shown in Figures [T] and |2] are snapshots of 
internal energy density and mass density at f = Jnr ~ 970 days. 
It should be noted that the GRB outflow at this moment is still 
highly anisotropic in the angular distribution of energy. 

Figure|3]shows the time evolution of the opening angles in- 
side which a certain percentage of the total energy, excluding 
rest mass energy, resides. At f = « 373 days, 50%, 90%, 
95% and 99% of the total energy are inside an opening angle 
of 0.14, 0.24, 0.28 and 0.41, respectively. At f = Tnr ~ 970 
days, 50%, 90%, 95% and 99% of the total energy is within 
an opening angle of 0.25, 0.54, 0.63 and 0.77, respectively. 
At f = 20 yeai-s, 50%, 90%, 95% and 99% of the total energy 
is within an opening angle of 0.72, 1.28, 1.43 and 1.55, re- 
spectively. Note that for a spherical blast wave, the opening 
angles would be 1.05, 1.47, 1.52 and 1.56, for 50%, 90%, 
95% and 99% of the total energy, respectively. It is clear that 
the spreading of energy to large angles happens very slowly. 
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Fig. 1. — Time evolution of the internal energy density in the local rest frame. The internal energy density (ergcm^^) is color coded on a logarithmic scale. 
Snapshots of the simulation are shown at (a) 147 days, (b) 256 days, (c) 372 days, (d) 970 days, (e) 27.6 years and (f) 150.3 years in the lab frame after the 
explosion. The beginning of the simulation is at 147 days. Panel (d) shows that the GRB outflow is still highly anisotropic axt = fa 970 days. The nearly 
vertical shock front near the equator in panels (e) and (f) is a Mach stem, which is a result of the shock collision along the equator The minimal value in the color 
scale corresponds to < 10"* ergcm"^. Note that the surrounding medium indeed has an internal energy density of ~ 2.25 X 10"'^ ergcm"^^. 
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Fig. 2. — Time evolution of the density in the local rest frame. The density (gem"') is color coded on a logarithmic scale. Snapshots of the simulation are 
shown at (a) 147 days, (b) 256 days, (c) 372 days, (d) 970 days, (e) 27.6 years and (f) 150.3 years in the lab frame after the explosion. The beginning of the 
simulation is at 147 days. Panel (d) shows that the GRB outflow is still highly anisotropic at f = ?nr 970 days. The vortex feature in panels (b) and (c) is an 
indication of the Kelvin-Helmholtz instability due to velocity shear The nearly vertical shock front near the equator in panels (e) and (f) is a Mach stem, which 
is a result of the shock collision along the equator. 



5 




1000 

Time (day) 



10000 



Fig. 3. — Time evolution of the jet opening angle. The opening angles 
inside which a certain percentage of the total energy, excluding rest mass 
energy, resides are shown for 50% (solid line), 90% (long dashed line), 95% 
(dotted line), and 99% (dashed line) of the total energy. Also shown in dash- 
dotted line is an analytic formula, 9 = 6oe^""'o)/'sNT 
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Fig. 4. — Angular distribution of energy. Different lines are for different 
times: t = (solid line), tg 373 days (clash-dotted line), ?nr ~ 970 days 
(dash dot dot line), 5fNR ~ 13 years (dotted line) and 150 years (dashed line). 



It takes 10-20 years for the jet to be somewhat spherical. 
The approach to a spherical blast wave and the Sedov-von 
Neumann-Taylor solution is a very slow process. 

The angular distribution of the total energy, excluding rest 
mass energy at f = 0, ~ 373 days, Tnr ~ 970 days, SJnr ~ 
13 years and 150 years is shown in Figure |4] The outflow 
evolves from a jet with an opening angle of = 0.2 into an 
isotropic blast wave. At f = 5fNR ~ 13 years, the energy per 
unit sold angle varies by about an order of magnitude between 
the axis and the equator During the evolution, the angular 
energy distribution is not in any kind of universal profile. In 
particular, the outflow does not resemble a top-hat jet. 

The very little sideways expansion we find is in agree- 
ment with the res ults of previous num erical simulations 
(iGranot et a l. 2001: Cannizzo et al.ir2004 . Figure |3] shows 
that the opening angle grows logarithrnically over time. How- 
ever, analytic work jRhoadsl 119991; ISari. Piran. & HalpernI 
Il999l) has predicted that the jet opening angle should grow ex- 



ponentially, 9 ^ eC('-'o)/'sNT (Fig.|3]l. Thus, according to these 
analytic estimates, the transition from jet-like to spherical-like 
takes place over practically no time. Why does the jet spread 
so rapidly in the analytic work compared with that in the nu- 
merical simulations? The main reason is that the jet in ana- 
lytic work is assumed to have an unrealistic top-hat distribu- 
tion as a function of angle. As it expands sideways, the top-hat 
jet has more and more working surface, which in turn rapidly 
decreases the Lorentz factor Since the speed of the sideways 
expansion depends upon the Lorentz factor, the sideways ex- 
pansion of a top-hat jet becomes a runaway process, which 
grows exponentially. In fact, the jet is far from top-hat dur- 
ing the sideways expansion (Fig. HJi. The information of the 
existence of surrounding medium outside the jet opening an- 
gle propagates towards the axis as a rarefaction wave, which 
moves at the sound speed in the local rest frame. The part of 
the jet that the rarefaction has not reached will behave exactly 
like a spherical outflow and continue to expand radially. The 
part of the jet that is affected laterally by the jet edge is in- 
side a rarefaction wave, which has a more complicated angu- 
lar profile than a top-hat. Moreover, the reverse shock which 
appears in the early stages further slows down the sideways 
expansion. Can numerical simulations underestimate the rate 
of sideways expansion? We believe that it is unlikely. In fact, 
numerical simulations tend to overestimate the rate of side- 
ways expansion of relativistical ly moving materia l due to the 
inevitable numerical viscosity dZhang & MacFadv en 2006). 
Also the Blandford-McKee solution is very challenging for 
numerical simulations due to its extremely thin structure be- 
hind the shock. The Lorentz factor of the material just be- 
hind the shock tends to be lower than the analytic Blandford- 
McKee solution (Fig. |5]l. This will also increase the rate of 
sideways expansion numerically. Hence, we conclude that the 
runaway lateral expansion derived in the analytic work does 
not exist in reality and the sideways expansion is a slow pro- 
cess. 

Figure |5] shows the time evolution of various properties of 
the fluid just behind the forward shock: the position of the 
shock front, the product of the Lorentz factor and radial ve- 
locity, and the internal energy density at various angles: 9 = 0, 
0. 19, and 7r/4. Also shown are the Blandford-McKee solution 
and Sedov-von Neumann-Taylor solution for comparison. It 
is shown that the outflow can be approximately described by 
the Blandford-McKee solution at early times (f < Jnr « 970 
days) and the Sedov-von Neumann-Taylor solution at late 
times (f > 5fNR ~ 5000 days). And the transition takes place 
over a period of ^ 4fNR. It is striking that the Blandford- 
McKee solution is valid for the material near the jet axis until 
the time when the Lorentz factor decreases almost to 1 . Note 
that the assumption of ultrarelativistic velocity upon which 
the Blandford-McKee solution is based has become invalid 
before that time. It should also be noted that even with 16 
levels of mesh refinement, the simulation still suffers from in- 
sufficient resolution at early times. 

The results of our hydrodynamic simulation are summa- 
rized as follows: (1) The initial condition of the jet is de- 
scribed by the Blandford-McKee solution; (2) The jet slows 
down as it sweeps up the surrounding medium; (3) At f ^ 
tg, the jet starts to undergo sideways expansion, but at a 
rate much slower than the exponential growth predicted by 
analytic work; (4) The jet becomes nonrelativistic and the 
Blandford-McKee solution breaks down at t ^ Tnr; (5) The 
outflow becomes more and more spherical and undergoes a 
slow transition into the Sedov-von Neumann-Taylor solution; 
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Fig. 5. — Properties of the fluid just beliind tlie forward shock as a function 
of time. The position of the forward shock, the product of the Lorentz factor 
and radial velocity, and internal energy density are plotted in Panels (a) {top), 
(b) {middle), and (c) {bottom), respectively. Properties at three different an- 
gles are shown: (1) S = {solid lines), (2) = 0.19 {dotted lines), and (3) 
6 = 7^1 A {dashed lines). In Panel (b), the unit for velocity is the speed of light. 
The Blandford-McKee solution and Sedov-von Neumann-Taylor solution are 
also shown in thick gray dashed lines. In Panel (b), instead of the product of 
the Lorentz factor and velocity, the Lorentz factor and velocity are plotted, 
for the Blandford-McKee solution and Sedov-von Neumann-Taylor solution, 
respectively. 

(6) After t ^ SJnr, the outflow is close to spherical and can 
be described by the Newtonian Sedov-von Neumann-Taylor 
solution. 

3. AFTERGLOW RADIATION OF GRBS 

We now calculate the afterglow radiation of GRBs using 
the data from our two-dimensional special relativistic hydro- 
dynamic simulation. 

3.1. Frames and Times 

Two frames are involved in the hydrodynamic simulation 
(§ 12.2b : the local rest frame of a fluid element and the lab 
frame of the GRB central engine. The measurements of 
spacetime in the above two frames satisfy the Lorentz trans- 
formation. For an observer on the earth, there is an additional 
frame: the observer frame ^ If a photon is emitted at time t 
and position r in the lab frame, an observer will receive it at 

?e = (l+z)(f-— ) (7) 
c 

after the main burst. Here, z is the redshift of the GRB, and 
n is the unit vector pointing from the burster to the earth. 
Clearly, photons received by the observer at a given time are 
emitted by different regions of the GRB outflow at different 

' The observer frame is th e same inertial frame as th e lab frame if cosmo- 
logical effects are neglected jZhang & Meszarosll2004l) . 




1000 10000 
Emission Time (day) 



Fig. 6. — Observer time vs. Emission time. The results for fluid elements 
just behind the shock at various angles are shown: 9 = {solid line), 0.19 
{dotted line), and 7r/4 {dashed line). The relations /gj = (l+z)//47^ {long 
dashed gray line) and = (1 +z)t/\(r)^ {dash dot dot gray line) are plotted 
for comparison. Here, 7 is assumed to obey the Blandford-McKee solution, 
and the cosmological redshift is set to z = 1. Also plotted is (9 = (1 +z)t 
{dash-dotted gray line). 

times. However, the emission is mostly from a very small re- 
gion if the outflow is relativistic. Therefore, a simple relation 
is often used in analytic work for the emission time and ob- 
server time. For a GRB outflow in the relativistic phase, the 
two times satisfy 

(8) 

here 7 is the Lorentz factor of the fluid just behind the shock, 
which is described by the Blandford-McKee solution. The 
factor of 4 in the denominator of Equation[8]is due to the effect 
that "typical" photons seen by an observer are emitted from 
off-axis regions (e.g.- .Waxm an 1997a). When the velocity of 
the outflow is Newtonian (i.e., v <C c), the relation is simply 

fe«(l+z)f. (9) 

There is, however, no good analytic relation for the transrela- 
tivistic phase. 

Our numerical approach has the advantage of treating the 
various times accurately. Figure |6] shows the observer time 
versus emission time for fluid elements just behind the shock 
at various angles. Analytic relations are also plotted for com- 
parison. We have discussed various timescales in § 12.21 in- 
cluding te, ts, tsm and Jnr (Eqs. [3] IH |5]and|6]l. All these 
timescales are measured in the lab frame. A good approx- 
imation is that Eq. [8] is valid for tg. The time consists 
of two parts (Eq. |4]i. Since the outflow is relativistic before 
tg and is then assumed to beco me Newtonian ve r y qui ckly 
due to the sideways expansion, iLivio & WaxmanI (120001) ar- 
gued that the GRB outflow will become spherical at fj ~ 

(1 +z)Rg/c - 372(1 +z)Ell%n~Q^\0o/Q.2f/^days. The time 

1/3 -1/3 

tNR.B ~ (1 +z)tNR ~ 970(1 +z)£iso,53"o ^ays is also often as- 
sumed to be the beginning of the Sedov-von Neumann-Taylor 
phase (e.g., Pirani.2005i) . But it is clearly shown in Figure |6] 
that the relation % = (1 +z)f cannot be used at times fj and Tnr- 
It would be a mistake to treat them as the observer times (after 
a cosmological redshift correction.) Another mistake is that 
the outflow at these times is neither spherical nor Sedov-von 
Neumann-Taylor as we have shown in § 12.21 The transition to 
the spherical Sedov-von Neumann- Taylor solution takes place 



7 



over a rather long period (~ SJnr) in the lab frame. But, it 
turns out that the two mistakes somewhat compensate each 
other At % = (1 +z)fNR , the observed radiation is emitted 
by nonrelativistic material, which is undergoing the transition 
into the spherical Sedov-von Neumann-Taylor phase. How- 
ever, (l+z)fsNT ~ 116(l+z)(£;/2 X 10-**ierg)i/3no'''^days is 
certainly inappropriate as the observer time for the beginning 
of either the nonrelativistic phase or the Sedov-von Neumann- 
Taylor phase. 

3.2. Standard Afterglow Model 

Our calculation of the radiation is based upon the standard 
external shock model for GRB afterglow s (Meszaros & Reei 
[T997t IWiiers etanflQQTt IWaxmanlll997 b a: Sari et al. 199S). 
We use the formalism introduced by Granotetal. (1999), 
which works nicely with data from numerical simulations. A 
large number of data dumps from the hydrodynamic simula- 
tion are stored. The numerical GRB outflow consists of many 
small fluid elements AV = iTrr^sindArAO at a discrete set 
of lab times. The velocity, number density and internal en- 
ergy density of each fluid element is known at each lab time. 
A fluid element is assumed to exist over a period of At de- 
pending upon the interval between dumps. The observed flux 
density at given by Eq.|7]over a period of Af^ = (1 +z)At 
due to radiation emitted by the fluid element is given by. 



^F(i^(j),t(s)-- 



l+z P'jv') 
^irdl^\l-(3.n)- 



■Ay, 



(10) 



here di is the luminosity distance, 7 is the Lorentz factor, (3 is 
the dimensionless velocity of the fluid element, and P'{v') is 
the emitted energy per unit volume per unit frequency per unit 
time measured in the fluid rest frame, where the frequency 
v' is also measured in the fluid rest frame. The frequencies 
measured in the fluid rest frame and observed on the earth are 
related by ^ 

v' = {\+z)l{l-p-n)v^. (11) 

The total observed flux density as a function of frequency 
and observer time, F(z^0,f®), can be calculated by combin- 
ing AF(i/0,%) over all the fluid elements at all discretized 
times. 

We consider synchrotron emission, but ignore synchrotron 
self-absorption and inverse Compton scattering, for the sake 
of simplicity. To compute synchrotron emission of a fluid el- 
ement with a number density n and an internal energy den- 
sity e, we use the standard approach of assuming that = e^e 
and es = ese, here e^ is the energy density of radiation emit- 
ting electrons, es is the energy density of the magnetic field, 
and te and es are two parameters. A power-law distribution 
is assumed for the radiation emitting electrons: A^(7e) ~ 7"^, 
where 7,, is the Lorentz factor of electrons, and p = 2.5 is a 
constant parameter. We follow the work of ISari et"an (Il998h 
for calculating synchrotron emission. For our purposes, it is 
sometimes more convenient to use the lab frame time. We 
have taken into account the effect of electron cooling on the 
spectral power of synchrotron emission. The critical value of 
the Lorentz factor of electrons is computed using 



Aaregt ' 



(12) 



^ There is a typo in Eq. 4 of lGranot et alj jl999l) . Tiiere should be a factor 
of 1 + z in the first argument of P'. 



here nig is the mass of an electron, aj is the Thompson cross 
section, 7 is the Lorentz factor of the fluid element, and t is 
the time in the lab frame. The spectral power of synchrotron 
emission is assumed to be a broken power-law with three seg- 
ments separated by i^c, the cooling frequency, and I'm, the typ- 
ical frequency of electrons with the minimal Lorentz factor. 
Assuming that the outflow obeys the Blandford-McKee solu- 
tion and the emission time and observer time are related by 
Eq. [8] the two break frequencies in the observer frame are 
given by-* 

-2 



:L8 X 10"(l+z)-'e 



1 -3/2^^-3/2 



'0 



147 days 



Hz 



= 5.3 X W"(l + zr'^€ 



l/2,-3/2p-l/2 -K-1/2, 



^F~''^ w"'/"'' H7 
'^iso,53"0 



and 



= 6.7 X mi+z)-'efelEl,,n-^ y^^^ ^^^^ 
= L8xlO-(l+z)'/^.feX^o^3^e';;^Hz, 



-3/2 



(13) 



Hz 



(14) 



here f®^ is the observer time in units of days. The peak flux 
density is given by 



= 9mi+z)efEi, 



1/2 j- 

>,53«o 



mJy. 



(15) 



Here, a factor of 0.88 was added to the expression for 
fiy.max (Eq. 11 in ISari et all 11998) to reflect the fact that 
the synchrotron emission at the peak frequency has a 
smoothly curved sh ape rather than a broken power-law shape 
(iGranot et alJl999h . The two break frequencies become equal 
to the critical frequency 



z/,,e =2.9 X lO\l+z)-'e-''\-'Erl,^n-'/^Hz 



at 



: 3.6 xl0^ef6yX£3 days, 



(16) 



(17) 



which corresponds to 

f,,e = 3.4 X 10^(1 +z)4eXo,53nodays. (18) 

In the calculation present in this paper, we assume that = 
0. 1 , eg = 0. 1 , and the GRB is located at z = 1 . At the beginning 
of the simulation, to = 147 days, the two break frequencies 
are fcO.e = 2.8 x lO'^Hz and t',„o.© = 1.1 x lO'^^Hz. Hence 
the GRB outflow is initially in the fast cooling regime. The 
transition from the fast cooling to slow cooling regime takes 
place at 360 days in the lab frame, which corresponds to 6.8 
days for the observer, and the transition frequency is f,_0 = 
4.6xlO"Hz. 

3.3. Results of Afterglow Calculation 

A large amount of data from the hydrodynamic simulation 
needs to be stored for the calculation of afterglow radiation 
observed on the earth. Without an adequate number of data 
dumps, the afterglow calculation will not be very accurate. 
However, an unnecessary amount of data would be generated 
if we simply dumped data equally spaced in time because of 
the power-law time dependence of the blast wave. For our 
simulation, 3000 data dumps equally spaced in logarithmic 
time have been made. To test whether 3000 dumps is suf- 
ficient. We have performed low-resolution runs with 5000 



3 Our Eq.[T3]for i/c is 16 times smaller than Eq. 11 in lSari et alj fT99l) 
due to a factor of 4 differ ence in the expression for 7c between our Ea ll2l and 
Eq. 6 in lSari et all fT993) . 
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Fig. 7. — Multi-frequency lightcurves. The afterglow radiation from the 
forward jet is shown in soUd lines, whereas the radiation from both the 
forward jet and counterjet is shown in dashed lines. Flux density for var- 
ious frequencies from radio to X-ray is plotted: lO' (red), lO'" (green), 
10" (blue), 10'2 (cyan), lO'^ (magenta), lO'** (yellow), lO'^ (purple), lO"" 
(aqua), and lO" Hz (black). P lotted in black d o tted lines for c omparison 
are lightcurves with slopes from lSari et al.l n999t):'Rh oadj <1999t) . It should 
be noted that these analytic lines are arbitrary in magnitude. An analytic 
lightcurve for frequency at 1 GHz in the nonrelativistic phase iFrail et aH 
[2000; Livio & Waxman 2000) is also plotted in red dash-dotted line. The 
vertical dotted lines are at 7.9 and 340 days. 




v(Hz) 

Fig. 8. — Spectra at various times. Thin black lines are results of our 
calculation for different times: fg) = 0.6 (solid line), 8 (dotted line), 30 
(dashed line), 100 (dash-dotted), 1000 (dash dot dot line), and 10000 days 
(long dash line). Thick gray lines are fitted using analytic formulae for syn- 
chrotron emission from nonthermal electrons with a power-law distribution. 
tSari etal.,1998,) . 



dumps. We found almost no difference in the results of the 
afterglow calculation between using all 5000 dumps and us- 
ing only 3000 of the 5000 dumps. Therefore we conclude that 
the frequency of dumping data is sufficient. 

In this calculation, we use a cosmology with Ho = 
71kms"' Mpc~\ riM = 0.27, and = 0.73. Using the hydro- 
dynamic data (§ I2.2| | and the standard afterglow model (§ 13. 2b . 
we have calculated multi-frequency lightcurves (Fig.lT). The 
calculated flux density covers a wide frequency range: from 
lO^Hz (radio) to lO'^Hz (X-ray), and a wide range of ob- 
server time: from 0.5 day after the burst to ^ 27 years. The 
spectra at various times are shown in Figure |8] 




1 10 100 1000 10000 

Time (day) 



Fig . 9 . — Temporal decay index as a function of time. Assuming ^t 
the temporal decay index a for various frequencies from radio to X-ray is 
plotted: 10' (red), lO'" (green), lO" (blue), 10'^ (cyan), 10'' (magenta), 
lO'"* (yellow), 10'^ (purple), 10" (aqua), and 10''' Hz (black). Note that the 
indices for the last four lines (from lO'"* to lO'^) are alm ost ide n tical. Also 
plotted for comparison are analytic results of Sari et alj 119991) : IFrail et aTl 
<200Q[) . The vertical dotted lines are at 7.9 and 340 days. 

3.3.1. Jet Break 

It is shown in Figure [T] that there is an achromatic break 
in the lightcurves at ~ 8 days. Using the analytic work of 
ISari et al. ( 1999), we can estimate that the jet break happens 

at tj ss 3.5(1 +z)£i!,/^53«o'''^(6lo/0.2)'*/3 days « 7.0 days. How- 
ever, we believe that while they happened to get the cor- 
rect answ er, t heir treatment of the sideways expansion was 
flawed (§ 12.2b . Furthermore, there is ambiguity in the rela- 
tion between the emission time and the observer time. The 
jet br eak would still exist even if there is no sideways expan- 
sion ("Meszaros & Reesl fl999l: iPanaitescu & Meszaroslll999l: 
lOranot & Kumar 2003^. Once the observer sees the edge of 
the jet, the observed flux will decrease rapidly due to the 
obvious difference between a spherical blast wave and a jet 
and the sideways expansion will help decrease the flux. It 
is usually thought that the break due to missing flux would 
be shallower than the break due to sideways expansion. This 
is incorrect because th e ima ge of a GRB afterglow is limb- 
brighten (Granot et aLlll999h . We will discuss this in more 
detail later in this section. According to this line of argu- 
ment, we can also estimate the jet break time. For the sake 
of simplicity, one can think that the information about the 
existence of an edge propagates towards the axis as a rar- 
efaction wave, which moves at the local sound speed 
Therefore, the angle of the head of the rarefaction wave is 
&RF ^ &i)-l~^Cs/c ^ 6q-j~^/^/3. An observer on the axis 
will see the rarefaction when 7^1/ 0rf ^ 6*0 ' ( 1 + 1 / V3). The 
jet break time is around the moment when that happens. This 
leads to 

tj « 3.9(1 +z)Eli%n-'^' (^^ j ' days, (19) 

where /iiso.53 is the isotropic equivalent energy in units of 
lO^^erg. 

The analytic results of ISari et al.1 (119981 119991) : iRhoadsl 
d 19991) are also plotted in Figure |7] for comparison. For ex- 
ample, they predicted that the flux density at high frequencies 

(i^0 > f,,!,©) evolves as 

cx a'. (20) 
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Note that in our model the transition from the fast cooUng to 
the slow cooling regime takes place at ^ 7 days (§ 13.21 Fig.[8]l, 
which is roughly the same time as the jet break. For i/g < 10'^ 

Hz, the flux density evolves as ~ f^' before the jet break, de- 

yiating from the analytic result that it should evolve as ~ t]^^ 
dSari et al.l ll998') probably due to still insufficient resolution 
for ultrarelativistic material, and then becomes essentially flat 
until the frequency is above the typical frequency i^m,®- Then 
the flux density is slightly steeper than . For i', ^ ~ 10'^ 
Hz, the flux density is essentially flat before the jet break, and 
is steeper than afterwards. For vq^ > 10'^, the flux den- 
sity is flatter than t^''^'^'^^^^ before the jet break, probably also 
due to still insufficient resolution for ultrarelativistic material 
(7 > 10). Then it becomes much steeper than t'^'. 

The most notable difference between our simulation results 
and previous analytic results (ISari et al.ll 19991 : lRhoadsiri999l) 
is the sharpness of the jet break in our results. At radio fre- 
quencies, the post-break lightcurves after the spectral break 
can be fitted approximately by ^ t^^'. But at optical and X- 
ray frequencies, the post-break lightcurves can no longer be 
fitted by a simple power-law. If the power-law form f"" is 
used for the fit, the lightcurves have a varying temporal de- 
cay index as shown in Figure |9] The temporal indices for 
= 10''* Hz and above are almost identical as expected for 
frequencies above both the cooling frequency I'cS} and typi- 
cal fre quency t^m.m- N ote that a statistical study of pre-Swift 
bursts (IZeh et al.l2006 ) found that the post-break temporal de- 
cay index was in the range of 1.30-3.03. It is shown in Fig- 
ure |9] that the post-break temporal decay index can increase 
to 3 and then decreases to below the electron power-law 
index p. This is consistent with the argument that jet break 
is caused by the missing flux outside the jet opening angle. 
It should be noted that photons received by an observer at a 
given observer time are not emitted at the same time in the lab 
frame because of the difference in light travel time for differ- 
ent parts of the outflow. It takes longer for photons emitted 
from off-axis parts of the jet to reach the observer than pho- 
tons emitted on the axis at the same lab frame time. Thus, at 
a given observer time, the radiation from the off-axis part of 
the jet was emitted earlier when the jet energy density (which 
is decreasing as the jet decelerates) was higher and is brighter 
than that from the axis. In ot her words, the imag e of a GRB 
afterglow is limb-brightened (i Granot et al.lll999l) . Once the 
edge of the jet is seen by the observer, the flux density will 
decrease rapidly. Since the missing flux is brighter, the tem- 
poral decay i ndex will overshoot after the jet break (see also 
iGranoj l2007h . Because the afterglow image is more limb- 
brighte ned for frequencies above the cooling frequency than 
below (iGranot et al.lll999l) . the overshooting is more promi- 
nent at higher frequencies than lower frequencies. 

The spectra at various times are shown in Figure [8] They 
can be fit by broken power-law curves as expected. This con- 
firms the analytic analysis in § 13.21 The outflow is initially in 
the fast cooling regime. The two break frequencies become 
equal at ^ 8 days. Then the slow cooling regime is entered. 
The typical frequency v^^q continues to decrease over time. 
However the cooling frequency I'c,® increases slowly from 
- 10'2 Hz at 8 days to - 10'"* Hz at 10,000 days. 

3.3.2. Nonrelativistic Regime 

We have discussed in § l2.2l that Tnr in the lab frame marks 
the beginning of nonrelativistic regime. Using the relation 
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Fig. 10. — Estimated energy as a function of observer time. With observed 
flux density at 1/^ = 1 GHz, one can es tim ate the total energy of the GRB 
outflow using Ea. \23\ (solid line) or Eq. |22] The latter estimate requires an 
estimate of the Sedov-von Neumann-Taylor time Tsnt- Thi'ee different values 
for ?sNT are used: 116 (dotted line), 232 (dashed line), and 58 da.ys(dash- 
dotted line). The true energy is 2 X 10^' erg. 

^ {\ +z)t/4-j^, we can obtain 

fNR,e « 170(1 +z)£,l£3«o'/'days. (21) 

When the GRB outflow becomes a Sedov-von Neumann- 
Taylor blast wave at lat e times, the hghtc u rve will become 
flatter (iDai & Lul [T999I: iFrail et al.l l2000t iHuang & Chenil 
l2003h . At late times, the radio frequency is typically in 
between 1/^.® and fc,© (Fig. [8]l. The radio lig htcurve in 
the Sedov-von Neumann-Tay lor phase evolves as (iFrail et al.l 
l2000HLivio & Waxmanll200a) 

F ~n9('l4-,^(3-p)/2/ \ / V^'^J/'^p j-2 ,-(/'-1)/2 
I'^KO.Iil+zy Vol/ vol/ " ^i.5ldL.28^^GHz 

/ X (21-I5p)/10 

here ii/, 28 is the luminosity distance in units of 10^** cm, and 
t'GHz is frequency in units of GHz. 

It should be emphasized that the GRB outflow at fNR can- 
not be described by the Sedov-von Neumann-Taylor solution 
yet because it is still highly nonspherical. In fact, it does not 
become approximately spherical until ^ SJnr- Since the tran- 
sition from relativistic to nonrelativistic flow is very smooth 
(§ I2.2I 1. it is reasonable that the lightcurves at late times be- 
come flat gradually (Figs.|7]and|9]l. There are no sharp breaks 
in them except the bumps from the c o unter jet at very lat e 
times. According to ' Frail et al.1 (l2000h : iBerger et all (|2004 : 
|Fiailetal.(2005) since the hghtcurve can be observed months 
or even years after the burst when the outflow is putatively 
a spherical Sedov-von Neumann-Taylor blast wave, the radio 
afterglow would allow for an accurate estimate of the total en- 
ergy of the GRB outflow. They used the flattening present in 
the observed radio afterglow lightcurves to estimate fsNX, and 
performed calorimetric analysis of the explosions. However, 
there is no sharp transition in the lightcurves calculated from 
our simulation (Fig. |7]i. It would therefore be very inaccurate 
to determine tsm from the rather smooth late-time afterglow. 
It is shown in Figure[TO]that the estimated energy using Eq.l22] 
could have an order of magnitude difference if the estimated 
fsNT varies by a factor of 4. Here, other parameters such as 
eg, eb, p and no are assumed to be known. The emission from 
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the counter jet can be observed at late times too. This causes 
a bump in the lightcurves (Fig.|7|. Eq.l22land our calculation 
do not match very well until ^ 50(l+z)fsNT ~ lO"* days. 
Fortunately, by substituting Eq.|5]into Eq.|22l we can obtain 
an equation without ts,m, 

f.«0.16(l+z)(^-">/2(^) (iL)'^\P9-.0,)/20^(5,.3)/10 
/ X (21-15rt/I0 

(9^) "^J^- ^''^ 

Here, a factor of 0.8 is added to the expression to fit the late- 
time afterglow better Figure [TO] shows that the estimated en- 
ergy using Eq.|23]is more accurate. Even at f ( 1 +z)fsNT /2 ~ 
100 days, it only overestimates the energy by ~ 50%. And the 
inferred energy using Eq. |23] agrees extremely well with the 
true energy after 10"* days. 

Fig. I2] shows that there is a bump in the lightcurves due 
to radiation from the counter jet. The late-time bump from 
the counter jet ha s been discussed by Granot & Loeb (2003); 
iLi & Sond (120041) : IWang. Huang. & Kong (2009). However, 
the bu mp in our results is rn u ch smoother and broader than 
that in lGranot & Loebl ([2003'); 'Li & Song" ('2004'). Our results 
are consistent with those of Wang et al. (2009). Initially, the 
counter jet moves relativistically. Then it slows down and be- 
comes nonrelativistic at Jnr in the lab frame. Its radiation can- 
not be observed by the observer on the earth until it becomes 
nonrelativistic at Jnr in the lab frame. This would lead to an 
estimate of the time for the counter jet bump in the lightcurve 

fcj,® ~ 2(1 +z)fNR ~ 1900(1 +z)Ell%n-''' days, (24) 
and an estimate of the ratio of the two fluxes at the peak 

P /, X (21-15p)/10 

for p = 2.5. In our calculation of the afterglow lightcurves, 
for = I GHz, the flux from the counter jet becomes com- 
parable to that from the forward jet at ^ 1700 days. At 3800 
days, the ratio of the two fluxes at 1 GHz reaches its peak of 
6. The above estimates (Eqs.l24landl25Tl agree extremely well 
with the numerical results. Because Jnr is the nonrelativistic 
timescale in the lab frame for the fluid on the jet axis, it is ex- 
pected that the first appearance of the counter jet bump from 
off-axis emission is earlier than fcj 0. And numerical calcula- 
tion indicates that the counter jet bump starts at ~ (1 + z)fNR. 

4. DISCUSSION 

Our relativistic h ydrodynamic simulation shows, as did 
iGranot et al.l (1200 lb . that the sideways expansion of a rela- 
tivisti c GRB jet is a very slow process (§ 12. 2t . Analytic works 
(e.g., Rhoadslll999l;ISari et al.lll999.) based on a top-hat distri- 
bution of energy greatly overestimated the rate of the side- 
ways expansion. A frequently asked question is, what is the 
speed of the sideways expansion? We think the question itself 
is wrong because it implicitly assumes a top-hat distribution 
during the expansion. This would always lead to an expo- 
nential growth of the jet opening angle. The jet during the 
expansion is indeed far from a top-hat distribution (Fig. |4|l, 
and it expands slowly. Another question is, what causes the 
jet break? Our calculations show that the jet break in GRB 
afterglow lightcurves is mainly caused by the missing flux 
when the edge of the jet is observed (§ 13.3. lb . Fortunately, the 



widely used formula, Eq. 1 in lSari et"al] (Il999h . which relates 
the jet break time to the properties of GRB jets is accurate. 

It is generally believed that the spherical Sedov-von 
Neuma nn-Taylor self-similar solution can be applied at ^ Jnr 
(e.g., Piran' l2005h . We, however, find that the time Jnr is 
the beginning of a rather slow transition into the Sedov-von 
Neumann-Taylor solution, and the outflow at that time is still 
highly nonspherical (§ 12.2b . Our simulation shows that the 
outflow can be described by the Sedov-von Neumann-Taylor 
solution after ^ 5rNR- Note that the time Tnr measured in 
the lab frame is not equivalent to the observer time since the 
velocity of the flow is not sufficiently small yet compared 
with the speed of light. Fortunately again, the afterglow at 
f© ~ (1 +z)fNR is during the nonrelativistic phase. 

We have found that the lightcurve after the jet break will 
become increasingly flatter over the time. But the flattening is 
a very gradual process. There is no characteristic timescale at 
which the lightcurve will suddenly become flatter Our results 
disagree with the common notion that the flux density evolves 
as c>c t'^ and then switches to oc f^' '^''^/'^ after Jnr or f,. 

However, late-time radio observations can reveal a wealth 
of information about the GRB outflow. Eq. |23] can be po- 
tentially useful in determining the true energy of the out- 
flow. We predict a late-time bump in flux density due to the 
radi at ion from the counter jet part of the outflow (see als o 
iGranot & Loebl I200I iLi & Sond l2004t IWang et al.l l2009l) . 
The radio afterglow of GRB 030329 had been observ ed up 
to 1128 days after the burst han der Horst et alJl2008l) . No 
firm conclusion was drawn as to whether or not the emission 
from the counter jet had been observed. However, the counter 
jet bump may have been observed already in GRB 980703, 
which was rnonito red in the radio band up to ^ 1000 days 
dBerger et alj|200l]) . It is reasonable to attribute the late-time 
radio flux to the host galaxy. However, it is also possible that 
the flux at ^ 1000 days was actually emitted by the counter 
jet. This possibility could be easily tested through a future 
radio observation of the host galaxy of GRB 980703. 

At the end of the hydrodynamic simulation (f sa 150 years), 
the GRB outflow is almost a sphere with an aspect ratio of 
~ 0.8. Thus, it would be very difficult to distinguish such 
a GRB remnant at an age of more than ~ 200 years from a 
supernova r emnant using morphology alone. However, previ- 
ous st udies (I Aval & Piranll200ll : lRamirez-Ruiz & MacFadyenI 
I2OO8I) have found that a GRB remnant in a similar interstellar 
medium will not become a sphere until ~ 5000 years. The 
striking difference is due to different initial setups. In our 
hydrodynamic simulation, the GRB outflow initially moves 
along the radial direction of spherical coordinates, whereas it 
moves parallel to the symmetric axis in their simulations. Ob- 
viously, the outflow sweeps up the surrounding medium at a 
much higher rate in our simulation than in theirs. Therefore, 
it is not surprising that the GRB outflow becomes a sphere at 
much earlier time in our simulation than in theirs. 

High-resolution hydrodynamic simulations of GRB out- 
flows are very computationally expensive. An intermediate 
approach to model the hydrodynamic evolution of GRB out- 
flows could be as follows. A nonspreading jet obeying the 
Blandford-McKee solution can be assumed until Tnr- For 
t > 5fNR, the outflow can be assumed to obey the Sedov-von 
Neumann-Taylor solution. Interpolation can be used for the 
transition phase Tnr < t < 5fNR (see Fig.|5]). 

Many aspects of our numerical calculations can be im- 
proved. The calculation of radiation could include syn- 
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chrotron self-absorption, which can be important for radio af- 
terglow, and the inverse Com pton process, which can be im- 
porta nt for X-ray afterglow (S ari & Esinll2001l: [Harrison et alj 
I2OOI . In our numerical simulation, the hydrodynamic evo- 
lution of the GRB outflow is adiabatic. This is a good ap- 
proximation provided that Eg-^l during the early fast cooling 
regime. However, radiative loss could greatly affect the hy- 
drodynamics of the outflow during the early afterglow phase 
(e.g., less than a day) if is close to 1 jCohen. Piran. & Saril 
fl998h . 

A uniform jet described by the Blandford-McKee solution 
is assumed as the initial condition of our numerical sim- 
ulation. This is justified because the jet initially moves at 
ultrarelativistic speeds and sideways expansion at very early 
times is likely to be modest. Besides uniform jet models 
based on the Blandford-McKee solu t ion, structured jet 
models (Meszaros, Rees, & Wiiers 19981; iPerna. Sari. & Fraill 
I2OO3I: IGranot & Kumar. ,2003.) have also been proposed to 
expla in various phenomena, in cl uding je t break ( Dai & Go3 



2001 



iRossi. Lazzati. & ReesI 120021: IZhang & MeszarosI 
Moreover, numerical simulations have shown that 
nonuniform structures are ex pected for relativistic jets emerg- 
ing from massive stars ("Zhang. Woosl ey. & MacFadvenl 



20021) . 



2001 



Zhane. Woosley. & Heger 



20041 



Morsonv. Lazzati. & Beg elmani |2007D . Recently 'Gruzinov| 
P007) has shown that the spherical Blandford-McKee solu- 
tion is not an attractor for a generic asymmetric explosion. 
Therefore, it is very important to investigate non-Blandford- 
McKee models and their consequences on GRB afterglows 
with high-resolution numerical simulations, which we will 
present in future publications. We speculate that the shallow 
deca y phase in the early X-ray afterglows of Swift bursts 
(e.g.. iNousek et al.ll200a lO'Brien et al.ll2006) could be due 
to non-Blandford-McKee behavior of the outflow and the 
normal decay phase is reached as the outflow approaches the 
Blandford-McKee solution. 

Magnetic field is not included in the two-dimensional rela- 
tivistic hydrodynamic simulation presented in this paper, be- 
cause we are mainly concerned with the stage when the mag- 
netic field is no longer dynamically important. However, since 
GRB jets m ay be pow ered by electromagnetic p rocesses (e.g.. 



Blandford & Znaiek 197 71 lUzden sky & MacFadyen 2006 
2007; Barkov & Komissarovl 120081: .Bucciantini et al.. .200S| 



magnetic field might have important effects on the dynam- 
ics of the jet during its pre- Blandford-McKee phase and 
early afterglows (e.g., IZhang & Kobayashi 2005). Prelimi- 
nary one-dimensional relativistic magnetohydrodynamic sim- 
ulations have found that the early afterglows are strong depen- 
dent o n the magnetization of the GRB outflow dMimica et alj 
l2009h . 

The environ ment of a GRB can have a huge impact on its 
afterglow (e.g.. Me szaros etal.lll998t IChevalier & OI1999L 
I2OOO.) . Thus, another important issue is to investigate GRB 
outflows propagating in stellar winds since the pr ogenitors 
of long soft GRBs are believe d to be ma ssive stars (" Woosleyl 
[1993: Paczvnski 1998: MacFadven & Woosley 1999) . More 



recently, late-time Chandra observations of the X-ray after- 
glow of GRB 60729 up to 642 d ays after the burst have 
been reported bv lGrupe et alj (|2009|) . It is interesting that the 
lightcurve at such late times shows signs of steepening pos- 
sibly due to either jet break or spectral break rather than flat- 
tening due to nonrelativistic motion. This might indicate that 
the onset of the Sedov-von Neumann-Taylor phase has been 
decayed due to the wind medium environment 

We adopted the standard approach in modeling GRB after- 
glow radiation. For example, the magnetic field is assumed 
to carry constant fractions of the internal energy density, and 
so do synchrotron emitting nonthermal relativistic electrons. 
Moreover, the fractions are assumed to be constant through- 
out the entire afterglow phase. The parameterization is a re- 
sult of the lack of understanding of these processes. It re- 
mains unclear how magnetic fields are generated and how 
particles are accelerated in GRB outflows. The Weibel insta- 
bihty is a plausible mech anism dGruzinov & Waxmanlll999l : 
iMedvedev & Loeblll999t) . Recent plasma simulations of the 
Weibel instability in relativistic collisionless shocks have 
shown promising results (Nish ikawa et al. 2003; Spi tkovskyl 
I2OO8I) . More recently, Zhan g. MacFadyen. & Wand (120091) 
have demonstrated amplification of magnetic field by a tur- 
bulent dynamo triggered by the Kelvin-Helmholtz instabil- 
ity with three-dimensional relativistic magnetohydrodynamic 
simulations. For con ditions relevant for late afterglow and 
prompt GRB emission lZhang et"an (l2009h obtained ^ 5 x 

io-\ 

In this paper, we have calculated afterglows for an ob- 
server on the axis of the GRB jet. However, orphan after- 
glows are expected for an off-axis observer, who has missed 
the main burs t because of the relativistic beaming effect 
(lRhoadslll997h . The rate of orphan af t erglows have been 
estimated analytically dOalal et a 1120021; iGranot et all 120021: 
i Levinson et al.ll2002l: " lNakar et al.ll2002lK Thus far, no orphan 
afterglow has been detected. A recent survey of 68 local 
Type Ibc supernovae, including 6 broad-lines supernovae also 
known as "hypernovae " found no evidence of off-axis GRBs 
dSoderberg et al.l2006h . Note that we have found that the side- 
ways expansion has been overestimated in previous analytic 
works. Thus, previous analytic works tend to overestimate 
the rate of orphan afterglows. The issue of off-axis GRB af- 
terglows will be investigated in a future publication. 
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